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METASTABILITY IN A STOCHASTIC NEURAL NETWORK 
MODELED AS A VELOCITY JUMP MARKOV PROCESS* 

PAUL C. BRESSLOFFt AND JAY M. NEWBY* 

Abstract. 

One of the major challenges in neuroscience is to determine how noise that is present at the 
r*f\ molecular and cellular levels affects dynamics and information processing at the macroscopic level 

of synaptically coupled neuronal populations. Often noise is incorprated into deterministic network 
models using extrinsic noise sources. An alternative approach is to assume that noise arises intrinsi- 
cally as a collective population effect, which has led to a master equation formulation of stochastic 
neural networks. In this paper we extend the master equation formulation by introducing a stochastic 
model of neural population dynamics in the form of a velocity jump Markov process. The latter has 
t*H the advantage of keeping track of synaptic processing as well as spiking activity, and reduces to the 

■^^ neural master equation in a particular limit. The population synaptic variables evolve according to 

piecewise deterministic dynamics, which depends on population spiking activity. The latter is char- 
Vj acterised by a set of discrete stochastic variables evolving according to a jump Markov process, with 

transition rates that depend on the synaptic variables. We consider the particular problem of rare 
transitions between metastable states of a network operating in a bistable regime in the deterministic 
limit. Assuming that the synaptic dynamics is much slower than the transitions between discrete 
spiking states, we use a WKB approximation and singular perturbation theory to determine the 
mean first passage time to cross the separatrix between the two metastable states. Such an analysis 
can also be applied to other velocity jump Markov processes, including stochastic voltage-gated ion 
channels and stochastic gene networks. 
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1. Introduction. Noise has recently emerged as a key component of many bi- 
ological systems including the brain. Stochasticity arises at multiple levels of brain 
function, ranging from molecular processes such as gene expression and the opening 
of ion channel proteins to complex networks of noisy spiking neurons that generate 
behaviour [28]. For example, the spike trains of individual cortical neurons in vivo 
tend to be very noisy, having interspike interval (ISI) distributions that are close to 
Poisson [72]. At the network level, noise appears to be present during perceptual de- 
cision making [79] and bistable perception, the latter being exemplified by perceptual 
switching during binocular rivalry [51, 71, 81]. Noise also contributes to the generation 
of spontaneous activity during resting states [23, 22]. At the level of large-scale neural 
systems, as measured with functional MRI (fMRI) imaging, this ongoing spontaneous 
activity reflects the organization of a scries of highly coherent functional networks 
• • that may play an important role in cognition. One of the major challenges in neuro- 

science is to develop our understanding of how noise that is present at the molecular 
and cellular levels affects dynamics and information processing at the macroscopic 
level of synaptically coupled neuronal populations. Mathematical and computational 
modeling arc playing an increasing role in developing such an understanding [45]. 

Following studies of biochemical and gene networks [74, 28], it is useful to make 
the distinction between intrinsic and extrinsic noise sources. Extrinsic noise refers 
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to external sources of randomness associated with environmental factors, and is of- 
ten modeled as a continuous Markov process based on Langevin equations. On the 
other hand, intrinsic noise typically refers to random fluctuations arising from the 
discrete and probabilistic nature of chemical reactions at the molecular level, which 
are particularly significant when the number of reacting molecules N is small. Under 
such circumstances, the traditional approach to modeling chemical reactions based 
on the law of mass action is inappropriate. Instead, a master equation formulation 
is necessary in order to describe the underlying jump Markov process. In the case of 
single cortical neurons, the main source of extrinsic noise arises from synaptic inputs. 
That is, cortical neurons are being constantly bombarded by thousands of synaptic 
currents, many of which are not correlated with a meaningful input and can thus be 
treated as background synaptic noise. The main source of intrinsic fluctuations at the 
single cell level is channel noise, which arises from the variability in the opening and 
closing of a finite number of ion channels. The resulting conductance-based model of a 
neuron can be formulated as a stochastic hybrid system, in which a piecewise smooth 
deterministic dynamics describing the time evolution of the membrane potential is 
coupled to a jump Markov process describing channel dynamics [63, 15, 58]. 

It is not straightforward to determine how noise at the single cell level trans- 
lates into noise at the population or network level. A number of methods involve 
carrying out some form of dimension reduction of a network of synaptically-coupled 
spiking neurons. These include population density methods [59, 61, 46], mean field 
theories [1, 35, 14, 13, 2], and Boltzmann-like kinetic theories [20, 67, 19]. However, 
such methods tend to consider either fully-connected or sparsely connected networks 
and simplified models of spiking neurons such as the integrate-and-fire (IF) model. 
Nevertheless, one interesting result that emerges from the mean-field analysis of IF 
networks is that, under certain conditions, even though individual neurons exhibit 
Poisson-like statistics, the neurons fire asynchronously so that the total population 
activity evolves according to a mean-field rate equation with a characteristic acti- 
vation or gain function [1, 35, 14, 13]. Formally speaking, the asynchronous state 
only exists in the thermodynamic limit N — > oo, where N determines the size of the 
population. This then suggests a possible source of intrinsic noise at the network 
level arises from fluctuations about the asynchronous state due to finite size effects 
[50, 48, 73, 6]; this is distinct from intrinsic noise at the single cell level due to chan- 
nel fluctuations and it is assumed that the latter is negligible at the population level. 
The presence of finite-size effects has motivated developing a closer analogy between 
intrinsic noise in biochemical and neural networks [7, 8], based on a rescaled version 
of the neural master equation introduced by Buice et. al. [17, 18], see also [60]. 

In the Buice et al master equation [17, 18], neurons are partitioned into a set 
of M local homogeneous populations. The state of the ath population at time t is 
specified by the number N a (t) of active (spiking) neurons in an infinite background 
sea of inactive neurons. (This is reasonable if the networks are in low activity states). 
Transitions between the states are given by a one-step jump Markov process, with the 
transition rates chosen so that standard Wilson-Cowan or activity-based equations are 
obtained in the mean-field limit, where statistical correlations can be ignored. One of 
the features of the Buice et. al. master equation is that there does not exist a natural 
small parameter, so that it is not possible to carry out a diffusion-like approximation 
using, for example, a system-size expansion. Indeed, the network tends to operate in a 
regime close to Poisson-like statistics. Neverthlcss, it is possible to solve the moment 
hierarchy problem using either path-integral methods or factorial moments [17, 18]. 
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In contrast, the Brcssloff master equation [7, 8] assumes that there is a finite number 
N of neurons in each local population and characterizes the state of each population 
in terms of the fraction of neurons N a (t)/N that have spiked in an interval of width 
At. The transition rates are rescaled so that in the thermodynamic limit N — > oo, 
one recovers the Wilson-Cowan mean-field equations. For large, but finite N, the 
network operates in a Gaussian-like regime that can be described in terms of an 
effective neural Langevin equation [7, 8] . One of the advantages of this version of the 
master equation from a mathematical perspective, is that a variety of well-established 
methods from the analysis of chemical master equations can be generalized to the 
neural case. For example, a rigorous analysis of the Langevin approximation can be 
carried out [16] by extending the work of Kurtz [43] on chemical master equations. 
Moreover WKB methods can be used to analyze rare transitions between metastable 
states, for which the Langevin approximation breaks down [8] . For a discussion of the 
differences between the two master equations from the perspective of corrections to 
mean field theory, see [76]. 

In this paper we go beyond the neural master equations by formulating the net- 
work population dynamics in terms of a stochastic hybrid system described by a "ve- 
locity" jump Markov process. This generalization is motivated by a major limitation 
of the neural master equations. That is, they neglect synaptic dynamics completely, 
only keeping track of changes in spiking activity. This implies, for example, that 
the relaxation time r for synaptic dynamics is much smaller than the fundamental 
time step At for jumps in the number of active neurons. Our model associates with 
each population two stochastic variables U a (t) and N a (t). The synaptic variables 
U a (t) evolve according to piecewise-deterministic dynamics describing, at the popu- 
lation level, synapses driven by spiking activity. These equations are only valid be- 
tween jumps in spiking activity, which are described by a jump Markov process whose 
transition rates depend on the synaptic variables. Formally speaking, the resulting 
stochastic dynamics can be modeled in terms of a differential Chapman-Kolmogorov 
(CK) equation: 

x -(v a (u,n)p(u,n,i))H V A(n,m;u)p(u,m,£). (1.1) 
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Here n = (n\, . . . ,um), u = (u\, . . . ,Um), and p(u, n, t) is the state probability 
density at time t. The drift "velocities" v a (u, n) for fixed n represent the piecewise- 
deterministic synaptic dynamics according to 

T^=v a (u,n), a = l,...,M, (1.2) 

and A represents the u-dependent transition matrix for the jump Markov process 
with J2 n ^( n i m i u ) — f° r au m - Note that the transition rates are scaled by a 
second time constant r a that characterizes the relaxation rate of population activity. 
In the limit r — >• for fixed r , equation (1.1) reduces to the neural master equation 
[17, 18, 7] with u = u(n) such that v Q (u(n), n) = 0. On the other hand, if r a — »■ for 
fixed r, then we obtain deterministic voltage or current-based mean-field equations 

T—jjr =v a (u, (n)), (1.3) 

where (n) = J^ n p(u, n)n with p(u, n) the steady-state density satisfying the equation 
Sm ^( n i m ; u )p( u j m ) = 0- It is straightforward to show using the Perron-Frobenius 
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Theorem that the steady-state density exists and is unique. Note that the limit r a — > 
is analogous to the slow synapse approximation used by Ermentrout [25] to reduce 
deterministic conductance-based neuron models to voltage-based rate models. Now 
suppose that the network operates in the regime < r a /r = e <C 1, for which there 
are typically a large number of transitions between different firing states n while the 
synaptic currents u hardly change at all. This suggests that the system rapidly con- 
verges to the (quasi) steady state p(u, n), which will then be perturbed as u slowly 
evolves. The resulting perturbations can be analyzed using a quasi-steady-state (QSS) 
diffusion or adiabatic approximation [64, 33, 57], in which the CK equation (1.1) is 
approximated by a Fokker-Planck equation. The latter captures the Gaussian-like 
fluctuations within the basin of attraction of a fixed point of the mean-field equa- 
tions, and can be used to investigate effects such as the noise-induced amplification 
of subthreshold oscillations (quasicycles) along smililar lines to [8]. However, the 
diffusion approximation for small e breaks down when considering rare event transi- 
tions between metastable states. (A similar problem arises in approximating chemical 
master equations by a Fokker-Planck equation in the large N limit [38, 24]). 

In this paper we will show how asymptotic methods recently developed to study 
metastability in stochastic ion channels and gene networks [58, 53, 54] can be extended 
to analyze metastability in stochastic neural networks. All of these systems are mod- 
eled in terms of a stochastic hybrid system evolving according to a CK equation of 
the form (1.1). For example, in the case of ion channels, n a would represent the 
number of open channels of type a, whereas u would be replaced by the membrane 
voltage v. The neural system is distinct, in that the numbers of discrete and contin- 
uous variables are equal. The structure of the paper is as follows. In §2, we present 
our stochastic network model and the associated neural CK equation, and carry out 
the QSS diffusion approximation for small e. We then analyze bistability in a one- 
population model (§3), and a two-population model consisting of a pair of excitatory 
and inhibitory networks (§4). In both cases, we carry out an eigenfunction expansion 
of the probability density and equate the principal eigenvalue with the inverse mean 
first pasage time from one metastable state to the other. The principal eigenvalue is 
expressed in terms of the inner product of a quasistationary density and an adjoint 
eigenfunction. The former is evaluated using a WKB approximation, whereas the 
latter is determined using singular perturbation theory, in order to match an absorb- 
ing boundary condition on the separatrix between the basins of attraction of the two 
metastable states. In the two-population case, calculating the effective potential of 
the quasistationary density requires identifying an appropriate Hamiltonian system, 
which turns out to be non-trivial, since the system does not satisfy detailed balance. 

A number of general comments are in order before proceeding further. 

(i) There does not currently exist a rigorous derivation of population rate-based mod- 
els starting from detailed biophysical models of individual neurons. Therefore, the 
construction of the stochastic rate-based model in §2 is heuristic in nature, in order 
to motivate the neural rate equations used in this paper. 

(ii) We use formal asymptotic methods rather than rigorous stochastic analysis to 
determine the transition rates between metastable states in §3 and §4, and validate 
our approach by comparing with Monte-Carlo simulations. Such methods have been 
applied extensively to Fokker-Planck equations and master equations as reviewed in 
[69], and provide useful insights into the underlying dynamical processes. In this 
paper we extend these methods to a stochastic hybrid system. One could develop a 
more rigorous approach using large deviation theory [30, 77], for example, although 



NEURAL VELOCITY JUMP MARKOV PROCESS 5 

as far as we are aware this has not been fully developed for stochastic hybrid sys- 
tems. Moreover, large deviation theory does not generate explicit expressions for the 
prcfactor, which we find can contribute significantly to the transition rates. 

(iii) We focus on networks operating in the bistable regime, where the simpler QSS 
diffusion approximation breaks down. There are a growing number of examples of 
bistability in systems neuroscience, including transitions between cortical up and 
down states during slow wave sleep [21, 65], working memory [37], and ambiguous 
perception as exemplified by binocular rivalry [3, 51, 44, 12]. On the other hand, 
in the case of oscillator networks, a diffusion approximation combined with Floquet 
theory might be sufficient to capture the effects of noise, including the noise-induced 
amplification of coherent oscillations or quasicycles [49, 4, 8]. An interesting issue 
is whether or not the WKB method and matched asymptotics can be applied to a 
network operating in an excitable regime, where there is a stable low activity resting 
state such that non-infinitesimal perturbations can induce a large excursion in phase 
space before returning to the resting state. One of the difficulties with excitable sys- 
tems is that there docs not exist a well-defined separatrix. Nertheless, it is possible to 
extend the asymptotic methods developed here to the excitable case, as we will show 
elsewhere within the context of spontaneous action potential generation in a model 
of an excitable conductance-based neuron with stochastic ion channels. 

2. Stochastic network model and the Chapman-Kolomogorov equation. 

Suppose that a network of synaptically coupled spiking neurons is partitioned into a 
set of M homogeneous populations with N neurons in each population, a = 1, . . . , M. 
(A straightforward generalization would be take to take each population to consist of 
0{N) neurons). Let \ denote the population function that maps the single neuron 
index i = 1, . . . , NM to the population index a to which neuron i belongs: x(*) = a - 
Furthermore, suppose the synaptic interactions between populations are the same 
for all neuron pairs. (Relaxing this assumption can lead to additional sources of 
stochasticity as explored in Ref. [29, 75]). Denote the sequence of firing times of the 
jth neuron by {T™, m € Z}. The net synaptic current into postsynaptic neuron i due 
to stimulation by the spike train from presynaptic neuron j, with \{i) = a > x(j) = Pi 
is taken to have the general form iV^ 1 J2 m ^afsit — T™ 1 ), where N~ 1 $ a p(t) represents 
the temporal filtering effects of synaptic and dendritic processing of inputs from any 
neuron of population f3 to any neuron of population a. For concreteness, we will take 
exponential synapses so that 

«M*) = w af) $(t), *(*) - T-h-^ T H{t) (2.1) 

(A more general discussion of different choices for $ Q /?(t) can be found in the reviews 
of Ref. [26, 9]). Assuming that all synaptic inputs sum linearly, the total synaptic 
input to the soma of the ith neuron, which we denote by Ui (t) , is 

u *w = E^ E *a/j(*-ir)=r E < m*-*')^ E a ^>' 

j;xU)=P °° /3 r,xU)=P 

(2.2) 

for all x(i) = a, where 

a,w = E^- T .n- ( 2 - 3 ) 
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That is, a,j(t) represents the output spike train of the jth neuron in terms of a sum of 
Dirac delta functions. (Note that in (2.2) we are neglecting any z-dependent transients 
arising from initial conditions, since these decay exponentially for any biophysically 
based model of the kernels $ Q ^). In order to obtain a closed set of equations, we have 
to determine threshold conditions for the firing times T™. These take the form 

rr = inf{M > Tr'mt) = « th , ^ > o>, (2.4) 

where « t h is the firing threshold and Vi(t) is the somatic membrane potential. The 
latter is taken to evolve according to a conductance-based model 

C^T = -Jcon.iW, • . •) + «i, ( 2 - 5 ) 

at 

which is supplemented by additional equations for a set of ionic gating variables [27]. 
(The details of the conductance-based model will not be important for the subsequent 
analysis). Let a a (t) denote the output activity of the crth population: 

and rewrite equation (2.2) as 

u > (*) = I Yj $q ^ - *>/j(*')* / > X(i) = a. 

J — oo a 



Since the right-hand side is independent of i, it follows that Uj(i) = u a (t) for all 
X{i) = a with 

m ft 

/3=1 J-oo 

In the case of exponential synapses (2.1), equation (2.7) can be converted to the 
differential equation 



du a 



M 

T "^f = -Ua(t) + ^2w a pap(t). (2.8) 

M /3=1 



In general, equations (2.2)-(2.5) are very difficult to analyze. However, consider- 
able simplification can be obtained if the total synaptic current m (t) is slowly varying 
compared to the membrane potential dynamics given by equation (2.5). This would 
occur, for example, if each of the homogeneous subnetworks fired asynchronously 
[34] . One is then essentially reinterpreting the population activity variables u a (t) and 
a a (t) as mean fields of local populations. (Alternatively, a slowly varying synaptic 
current would occur if the synapses arc themselves sufficiently slow [25, 10]). These 
simplifying assumptions motivate replacing the output population activity by an in- 
stantaneous firing rate a a (t) = F(u a (t)) with F identified as the so-called population 
gain function. Equation (2.7) then forms the closed system of integral equations 



/ ]T$^(t-iVMW- 

J —oo a 



u a (t)= >$ a0 (t-t')F(u {t'))dt'. (2.9) 
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The basic idea is that if neurons in a local population arc firing asynchronously then 
the output activity a a is approximately constant, which means that the synaptic 
currents are also slowly varying functions of time. A nonlinear relationship between 
a a and constant input current u a can then be derived using population averaging in 
order to determine the gain function F. One then assumes that the same relationship 
a a = F{u a ) also holds for time-dependent synaptic currents, provided that the latter 
vary slowly with time. In certain cases F can be calculated explicitly [1, 35, 14, 13]. 
Typically, a simple model of a spiking neuron is used, such as the intcgrate-and- 
fire model [34], and the network topology is assumed to be cither fully connected or 
sparsely connected. It can then be shown that under certain conditions, even though 
individual neurons exhibit Poisson-like statistics, the neurons fire asynchronously so 
that the total population activity evolves according to a mean-field rate equation with 
a characteristic gain function F. In practice, however, it is sufficient to approximate 
the firing rate function by a sigmoid: 

F(u) = ^ ,, (2.10) 

where 7, n correspond to the gain and threshold respectively. 

One of the goals of this paper is to develop a generalization of the neural master 
equation [17, 18, 7] that incorporates synaptic dynamics. We proceed by taking the 
ouput activity of a local homogeneous population to be a discrete stochastic variable 
A a (t) rather than the instantaneous firing rate a a = F(u a ): 

Mt) = *§. (2.11) 

where N a (t) is the number of neurons in the ath population that fired in the time 
interval [t — At,t], and At is the width of a sliding window that counts spikes. The 
discrete stochastic variables N a (t) are taken to evolve according to a one-step jump 
Markov process: 

N a (t) -)• N a (t) ± 1 : transition rate Q ± {U a (t), N a (t)), (2.12) 

with the synaptic current U a (t) given by (for exponential synapses) 



TdU a (t) = 



M 



-U a (t) + }^ Wal3 Ap(t) 

0=1 



dt. (2.13) 



The transition rates are taken to be (cf. [7]) 



n r 



Cl + {u a ,n a ) -¥ to+(u a ) = F(u a ), £l-(u a ,n a ) ->• fi_(n a ) = — . (2.14) 

The resulting stochastic process defined by equations (2.13), (2.11), (2.12) and (2.14) 
is an example of a stochastic hybrid system based on a piecewise deterministic pro- 
cess. That is, the transition rate fl + depend on U a , with the latter itself coupled 
to the associated jump Markov according to equation (2.13), which is only defined 
between jumps, during which U a (t) evolves detcrministically. (Stochastic hybrid sys- 
tems also arise in applications to genetic networks [82, 54] and to excitable neuronal 
membranes [63, 15, 41]). It is important to note that the time constant r a cannot be 
identified directly with membrane or synaptic time constants. Instead, it determines 
the relaxation rate of a local population to the instantaneous firing rate. 
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2.1. Neural master equation. Previous studies of the neural jump Markov 
process have effectively taken the limit r — > in equation (2.13) so that the continuous 
variables U a (t) are eliminated by setting U a (t) = Xifl w ai3^a(t)- This then leads to 
a pure birth-death process for the discrete variables N a (t). That is, let P(n,t) = 
Prob[N(i) = n] denote the probability that the network of interacting populations 
has configuration n = {ri\, ri2, ■ ■ ■ , »m) at time t, t > 0, given some initial distribution 
P(n, 0) with < n a < N. The probability distribution then evolves according to the 
birth-death master equation [17, 18, 7] 

' //>!i '' '•' ^ [(T a - 1) (fi-(n)P(n,t)) + (T- 1 - 1) (fi+(n)P(n,t))] , (2.15) 



rfi ^ 



where 



n+(n) = — P I ]T w a pn a /NM I , fi"(n) - ^, (2.16) 



/3 



and T a is a translation operator: T^ 1 P(n) = P(n a ±) for any function P with n a ± 
denoting the configuration with n a replaced by n a ± 1. Equation (2.15) is supple- 
mented by the boundary conditions P(n, t) = if n a = N + 1 or n a = — 1 for some a. 
The birth-death master equation (2.15) can be analyzed by adapting various meth- 
ods from the analysis of chemical master equations including system-size expansions, 
WKB approximations, and path integral representations [18, 7, 8, 16]. First, suppose 
that we fix At — 1 so that we obtain the Bressloff version of the master equation. 
Taking the thermodynamic limit N — > oo then yields the deterministic activity-based 
mean-field equation 

r a d ^ = -A a {t)+F{Y J W a0 A a (t)). (2.17) 

at 

(For a detailed discussion of the differences between activity-based and voltage-based 
neural rate equations, see Refs. [27, 9]). For large but finite N, the master equation 
(2.15) can be approximated by a Fokker-Planck equation using a Kramers-Moyal 
or system-size expansion, so that the population activity A a evolves according to a 
Langevin equation [7] . A rigorous probabilistic treatment of the thermodynamic limit 
of the neural master equation has also been developed [16], extending previous work on 
chemical master equations [42] . Although the diffusion approximation can capture the 
stochastic dynamics of the neural population at finite times, it can break down in the 
limit t —$■ oo. For example, suppose that the deterministic system (2.17) has multiple 
stable fixed points. The diffusion approximation can then account for the effects of 
fluctuations well within the basin of attraction of a locally stable fixed point. However, 
there is now a small probability that there is a noise-induced transition to the basin of 
attraction of another fixed point. Since the probability of such a transition is usually 
of order e~ cN , c = 0(1), except close to the boundary of the basin of attraction, such 
a contribution cannot be analyzed accurately using standard Fokker-Planck methods 
[78] . These exponentially small transitions play a crucial role in allowing the network 
to approach the unique stationary state (if it exists) in the asymptotic limit t — > oo, 
and can be analyzed using a WKB approximation of the master equation together 
with matched asymptotics [8]. In other words, for a multistable neural system, the 
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limits t — > oo and N — > oo do not commute, as previously noted for chemical systems 
[38]. 

Now suppose that we take the limit N — > oo,Ai — *• such that TV At = 1. 
We then recover the neural master equation of Buicc et. al. [17, 18]. In this case 
there is no small parameter that allows us to construct a Langevin approximation to 
the master equation. Nevertheless, it is possible to determine the moment hierarchy 
of the master equation using path integral methods or factorial moments, based on 
the observation that the network operates in a Poisson-like regime. The role of the 
sliding window size At is crucial in understanding the difference between the two 
versions of the master equation. First, it should be emphasized that the stochastic 
models arc keeping track of changes in population spiking activity. If the network 
is operating close to an asynchronous state for large N, then one-step changes in 
population activity could occur relatively slowly so there is no need to take the limit 
At —J- 0. On the other hand, if population activity is characterized by a Poisson 
process then it is necessary to take the limit At — > in order to maintain a one-step 
process. However, given the existence of an arbitrarily small time-scale At, it is no 
longer clear that one is justified in ignoring synaptic dynamics by taking the limit 
r — > in equation (2.13). This observation motivates the approach taken in this 
paper, in which we incorporate synaptic dynamics into the neural master equation. 
In the following, we will assume that the network operates in the Poisson-like regime 
in the absence of synaptic dynamics. 

2.2. Neural Chapman-Kolmogorov equation. Let us now return to the full 
stochastic hybrid system. Denote the random state of the full model at time t by the 
vector X(t) = {(U a (t) , N a (t)); a = 1, . . . , M}. Introduce the corresponding probabil- 
ity density 

Prob{U a (t) e (u a ,u a + du,N a (t) — n a :a — 1, . . . , M } = p(u, n, t\u , n , 0)du, 

(2.18) 
with u = (ui, . . . ,um) and x = (u, n). It follows from equations (2.13), (2.11), (2.12) 
and (2.14) that the probability density evolves according to the Chapman-Kolmogorov 
equation 

dp | 1 ^ 9[" tt (x)?(x,t)] . . 

-di + T^ &T a (2J9) 

ex. 

= - 2 [( T « - !) (w-(n Q )p(x,t)) + (T- 1 - 1) («+(« a )p(x,t))] , 

a 

with 

u)+(u a ) = F(u a ), u>-(n a ) = n a , v a (x) = -u a + ^ w af3 n P- (2.20) 

/3 

We have taken the limit N ->• oo, At ->• with NAt = 1. Note that equation (2.19) 
can be expressed in the general form of equation (1.1). Thus, in the limit r — > we 
recover the master equation of Buicc et. al. [17, 18], whereas in the limit r a — > we 
obtain the mean-field equations 

t-£- = (v a }{u(t)) = J2 «a( u (*), n)p(u(i), n) 

n 

M 

= -u a (t) + y~l w a p y~] rtff p(u(t) , n) . (2.21) 

/3=1 n 
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It can be shown that p(u, n) is given by a compound Poisson process with rates F(u a ), 
consistent with the operating regime of the Buice et. al. master equation [17, 18]. 
Hence, in this limit, 

(np) = F(up), (2.22) 

and (2.21) reduces to the standard voltage or current-based activity equation. 

2.3. Quasi-steady-state (QSS) diffusion approximation. In this paper, we 
will consider the regime in which the transitions between different firing states are 
much faster than the synaptic dynamics. Hence, fixing the units of time by setting 
t = 1, we take r Q /r = e <C 1. Since e<C 1, there will typically be a large number 
of transitions between different firing states n while the synaptic currents u hardly 
change at all. This suggests that the system will rapidly converge to the steady-state 
p(u, n) (if it exists) given by equation (2.24). The full probability density will then be 
perturbed away from this steady-state density as u slowly evolves. However, if e <C 1 
then these perturbations will be small and the solution will tend to remain close to the 
steady state. The resulting perturbations can then be analyzed using a quasi-steady- 
statc (QSS) diffusion or adiabatic approximation, in which the CK equation reduces to 
a Fokker-Planck (FP) equation. This method was first developed from a probabilistic 
perspective by Papanicolaou [64] , see also [33] . It has subsequently been applied to a 
wide range of problems in biology, including cell movement [62, 39], traveling- wavelikc 
behavior in models of slow axonal transport [68, 31, 32], and molecular motor-based 
models of random intermittent search [55, 57, 56, 11]. 

Consider a Chapman Kolmogrov equation of the general form (see equation (1.1)) 

— =-^2 g—{v a (u,n)p(u,n,t)) + -^A(n,m;u)p(u,m,t), (2.23) 

q=1 a m 

with e <C 1. We assume that for fixed u, the tensor A(n, m; u) is equivalent to 
a transition matrix. That is, suppose we relabel the discrete states according to 
n ->• I, m ->• J with I, J = 0, 1, . . . , \ f° r X — {N + 1) M and set A(n, m; u) = A u (u), 
p(u, n, t) = pi(u,t). Then the x x X matrix A(u) with elements Ajj(u) is taken to be 
irreducible and to have a simple zero eigenvalue with corresponding left eigenvector 
1 whose components are all unity. In other words, ^2j Ajj(u) = for all J. The 
Perron-Frobenius Theorem then ensures that all other eigenvalues are negative and 
the continuous-time Markov process for fixed u, 

dpi 1 v-^ . . 

— = -^A IjPj (u,t), 

has a globally attracting steady-state pi(u) such that Pi(u, t) —> p/(u) as t — >• oo. 
Here p/(u) is the unique right eigenvector corresponding to the zero eigenvalue of 
A(u), that is, ^ 7 Ajj(u)pj(u) = 0. In terms of the original notation, we have 

^^(n,m;u) = 0, J^ A(n, m; u)p(u, m) = 0. (2.24) 

n m 

In the following it will be convenient to introduce the summation operator 

[l T o/](u)=^/(u,n) (2.25) 
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for any function /(n, u). 

The first step in the QSS reduction is to decompose the probability density as 



p(u, n, t) — C(u, t)p(u, n) + ew(u, n, t), 



(2.26) 



where p(u, n) is the steady-state density given by equation (2.24), and 1 T o p = C, 
1 T o w — 0. Applying the summation operator to both sides of equation (2.23) and 
using J2 n M n > m, u) = gives 



dC 
~dt 



-l T o 



m 







52~-(v a [Cp + ew]) 



du a 



Next substitute (2.26) into equation (2.23) to give 

dC dw A y-^ d . ,^ 

_ p + e _ = A o«,-^-(„ a [Cp. 

a— 1 



where 



[A o w](u, n, t) = y^ A(n, m; u)w(u, m, t). 



(2.27) 



Combining with equation (2.27) shows that 

M 
= A n in -\- n~\ n 

dt 







E^r( w «[^+ ew ]) 



,a=l 



M 



^— K[c P + fW ]). 



a=l 



Collecting terms of leading order in e yields 



M 



= Y,j^^ c p)-p^ 



du, 



M 



Esr(^) 



9m, 



(2.28) 



The Frcdholm Alternative Theorem [66] ensures that this equation has a unique so- 
lution for w subject to the constraint 1 T o w — 0; we formally denote this solution 

as 



M 







AUJ2 7 ^(v«C P )-(Alop)l? 



du, 



M 







Esr(^) 



du a 



(2.29) 



where A^ is the pseudoinverse operator. Substituting for w back into equation (2.27) 
and ignoring C(e 2 ) terms finally gives the Fokker-Planck equation 



dC 
~dt 



M „ M M „ , „„ 

a=l/3=l Q v p 



du a 



where 



K, (u) = ^2 v <x ( u > n )/°( u ' n ) 



(2.30) 



(2.31) 
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and 

D a , 3 = ^K(m, u) - V a (u)]A^(m, n; u)[^(u, n) - V (u)]p{u, n). (2.32) 

m.n 

Note that we have expressed D a p in a more symmetric form using the fact that 
^ m A(m, n, u) = 0. Finally, let us introduce the function z a (u, m), which satisfies 
the equation 

^z Q (u, n)A(n, m; u) = -[v a (u,m) - V a (u)]. (2.33) 

n 

Since ^ m p(u, m)[v Q (u, m) — Va(u)] = 0, it follows from the Fredholm alternative 
that 

z a {u, n) = -^2 A \ m > n ; u)K(u, m) - V a (u)} (2.34) 



and thus 



D al 3 = ^ z a {u,n)[vfj{u,n) - V fj (u)]p(u,n). (2.35) 



Hence, under the QSS approximation, the stochastic dynamics is characterized by 
Gaussian fluctuations about the mean-field equations (2.21). However, as in the 
case of the system-size expansion of the neural master equation (2.15) for large N 
and At = 1, approximating the Chapman-Kolmogorov equation (2.19) by a Fokkcr- 
Planck equation for small e breaks down when considering rare event transitions 
between metastable states. This particular issue has recently been addressed within 
the context of stochastic ion channels and Hodgkin-Huxlcy dynamics [41], as well as 
gene networks [54], using asymptotic methods developed in [58, 53]. In the following 
sections, we will extend such methods to the neural CK equation. 

3. Metastable states in a one-population model. In order to develop the 
basic analytical framework, consider the simple case of a single recurrent population 
(M = 1) evolving according to the CK equation 



(3.1) 



dp d[v(u,n)p{u,n,t)] 
dt du 

~ -[uj+(u)p(u, n — l,t) + u;_(n + l)p(u, n + 1, t) 

-(u+(u) +u!-(n))p(u,n,t)], 

with boundary condition p(u, — 1, t) =0, drift term 

v(u,n) = —u + wn, (3-2) 

and transition rates 

to+(u) = F(u), uj-(n) = n. (3-3) 

Following the general discussion in §2.3, we expect the finite-time behavior of the 
stochastic population to be characterized by small perturbations about the stable 
steady-state of the pure birth-death process 

dp 1 

— = -[oj + (u)p(u,n - l,t) + w_(n + l)p(u,n + l,t) 

-(u+(u) +u-(n))p(u,n,t)}, 
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with u treated as a constant over time-scales comparable to the relaxation time of 
the birth-death process. The equation for the steady-state distribution p(u, n) can 
be written as [33] 

= J(u, n + 1) — J(u, n), 

with J(u,n) the probability current, 

J(u, n) = u)-(n)p(u, n) — oj + (u)p(u, n — 1). 

Since cj_(0) = and p(u,—l) = 0, it follows that J(u, 0) = and J(u,n) = for all 
n > 0. Hence, 

p(u, n) = ^\p(u, n-l)= p(u, 0) f[ ^ (3.4) 

u>- (n) - LJ - u)-(m) 

\ I m=1 \ ; 

with p(u, 0) = 1 - Em=i P( u > m )- 

Substituting the explicit expressions (3.3) for the transition rates, we have 

p(u,n)=p(u,0)f[ F M =p{U} o ) iE^L, (3.5) 

- LJ - m n! 

m— 1 

It follows that 

p(u, 0) = „«,!,„ ; , = e- F ^. (3.6) 

so the steady-state density is given by a Poisson process, 

[F( M )]"e- f ("> 

p(«,n) = j . (3.7) 

n! 

The mean number of spikes is thus (n) = F(u), and the mean-held equation obtained 
in the e — > limit is 

du „, . d^ 

— = -u + wf m = — — . (3.8 

at du 

The sigmoid function F(u) given by (2.10) is a bounded, monotonically increasing 
function of u with F(u) — > F as u — >• oo and -F(u) — >• as u — >• — oo. Moreover, 
F'(u) = 7.Fo/[4cosh 2 (7(u — ft)/2)] so that F(u) has a maximum slope at u — n given 
by 7F0/4. It follows that the function — u + wF(u) only has one zero if wjFo < 4 
and this corresponds to a stable fixed point. On the other hand, if w"(Fq > 4 then, 
for a range of values of the threshold k, [k\,k,2\, there exists a pair of stable hxed 
points u± separated by an unstable fixed point u* (bistability). A stable/unstable 
pair vanishes via a saddle-node bifurcation at K = K\ and k = k 2 - This can also be 
seen graphically by plotting the potential function ^(u), whose minima and maxima 
correspond to stable and unstable fixed points of the mean-field equation. An example 
of the bistable case is shown in Fig. 3.1. 

The problem we wish to address is how to analyze the effects of fluctuations 
(for < e <C 1) on rare transitions between the metastable states u± of the under- 
lying mean-field equation. As highlighted in §2.3, it is not possible to use a QSS 
diffusion approximation, since this only captures finite-time fluctuations within the 
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Fig. 3.1. Bistable potential \& /or i/ie deterministic network satisfying u — —u + F(u) = 
—dty/du with F given by the sigmoid (2.10) for 7 = 4, k — 1.0, i*b = 2. There exist two 
stable fixed points u± separated by an unstable fixed point u*. As the threshold n is reduced 
the network switches to a monostable regime via a saddle-node bifurcation. 

basin of attraction of a given metastable state. Therefore, we will proceed using 
the asymptotic methods recently introduced to analyze spontaneous action poten- 
tials in conductance-based single neuron model with stochastic ion channels [41]. 
As a first step, it is convenient to introduce the vector-valued probability density 
p(u,t) = (po(u,t), . . . ,p n (u, t), ■ ■ ■ ) with p n (u, t) = p(u,n,t). Equation (2.19) can 
then be rewritten in the matrix form 1 

| = -^(V( u )p) + iAp. (3.9) 

For given u, the diagonal drift matrix V has non-zero entries 

V n ^ n (u) = v n (u) = v(u, n) = —u + to, (3.10) 

and the tridiagonal transition matrix A has entries 

Ai,n-1 = <*>+(u), An >n = -u + (u) - w-(ri), A n . n+ i = w_ (n + 1). (3.11) 

Since A is a transition matrix, its columns sum to zero, it has one zero eigenvalue, 
and all other eigenvalues are negative. In particular, 

1 T A = 0, Ap = (3.12) 

with 1 T = (1, 1, . . .) and p n — p(u, n), where p(u, n) is the steady state density (3.5). 
In matrix notation, the mean-field equation (3.8) can be written as 

^=v(«).p(«). (3.13) 

1 Recall that we have taken the limit N —> 00, At — > such that Af At = 1. Hence, we are dealing 
with infinite matrices. However, this does not cause any problems. Indeed, we could equally well 
proceed by assuming that N is finite and At = 1/N, performing all calculations, and then taking the 
limit N — > co . The advantage of working with infinite N is that the steady-state density p is given 
by a Poisson process and we don't have to worry about boundary conditions at n = N. 
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Finally, note that since v n (u) = —u + wn, it follows that if the initial synaptic current 
u > then u(t) > for all t > 0. 

3.1. Quasistationary approximation. Suppose that the neural population 
starts in the left-hand well of the potential function *&(u), see Fig. 3.1, at the stable 
low activity state w_. On short time scales the solution rapidly converges to a quasi 
stationary solution that is only distributed across the left well. However, on a longer 
time scale, probability slowly leaks into the right well until the full stationary solution 
is reached. In order to estimate the exponentially small transition rate from the left 
to right well, we place an absorbing boundary at the unstable fixed point w*. (The 
subsequent time to travel from u* to the high activity fixed point u + is insignificant, 
and can be neglected). Thus, the CK equation (3.9) is supplemented by the absorbing 
boundary conditions 

p n {u*,t) = 0, forn = 0,...,fc- 1, (3.14) 

where < k < oo is the number of firing states for which the drift v(u*,n) < 0. The 
initial condition is taken to be 

Pn(u, 0) = 5{u - u_)5 n , no . (3.15) 

Let T denote the (stochastic) first passage time for which the system first reaches u* , 
given that it started at u_ . The distribution of first passage times is related to the 
survival probability that the system hasn't yet reached w* : 

5(t) = V/ Pn (u,t)du. (3.16) 

n=0 J ° 

That is, Prob{£ > T} — S(t) and the first passage time density is 

n=0 J ° 

Substituting for dp n /dt using the CK equation (3.9) shows that 

/(*)=>/ du= > v(u*,n)p n (u*,t). (3.18) 

We have used the fact that 1 T A = and p(0, t) — 0. The first passage time density 
can thus be interpreted as the probability flux J(us,,t) at the absorbing boundary, 
since we have the conservation law 

00 ft f) i °° 

^— ^ = - — , J(u,t) = ^2v(u,n)p n (u,t). (3.19) 

n=0 n=0 

The probability flux at the absorbing boundary can be approximated using a 
spectral projection method [80, 40, 58, 41]. Consider an cigenfunction expansion of 
the form 

DC' 

3=0 
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where the eigenfunctions satisfy the equation 



s< v «4 



L4, j s—(V4, i )--A<t, 1 =\ i <t, 1 , (3.21) 



together with the boundary conditions 

{4>j)n(u*) = 0, forn = 0,...,fc-l. (3.22) 

If the absorbing boundary is replaced by a reflecting boundary then there is a single 
zero eigenvalue Ao, for which the resulting stationary density p s is the corresponding 
cigenfunction. On the other hand, when there is an absorbing boundary, the station- 
ary solution ceases to exist due to a nonzero probability flux through u*. Moreover, 
Ao is perturbed away from zero but is exponentially small compared to the remaining 
eigenvalues. In other words, Ao = C(e _L / £ ) for some L > and Xj = 0(1), j > 1. 
It follows that all other eigenmodes decay to zero much faster than the perturbed 
stationary density. Thus, at large times we have the quasistationary approximation 

p( U ,t)~C c- Ao >o(w). (3.23) 

Substituting such an approximation into equation (3.18) implies that 

/(i)~vK)-0oK)Coe- Aot , Ai*»l, (3.24) 

The next step of the spectral projection method is to define a set of eigenfunctions 
for the adjoint operator, which satisfy the equation 



d _ , 1 
= - V 

and the boundary conditions 



£%' = - V ^) - ^ A % = A ^< ( 3 - 25 ) 



(&)„(«*) = 0, n > k. (3.26) 

The two sets of eigenfunctions form a biorthogonal set with respect to the underlying 
inner product, which is taken to be 

(f,g)= f ' f T (u)g(u)du. (3.27) 

Jo 

Now consider the identity 

(0o,£*So> = Ao<0o,So>. (3.28) 

Suppose that the exact eigenfunction </> satisfying the absorbing boundary conditions 
can be approximated by a quasistationary solution <fi e for which L(p e = without any 
absorbing boundaries. Under such an approximation, integrating by parts the left- 
hand side of equation (3.28) picks up a boundary term so that 

<ffK)VK)6)("») ,„ 9cn 

A — . (6.ZU) 

(<P e ,€o) 

The calculation of the principal eigenvalue Ao thus reduces to the problem of deter- 
mining the quasistationary density 4> e and the exact adjoint eigenfunction £o using 
perturbation methods (see below). Once Ao has been evaluated, we can then identify 
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the mean first passage time (T) with A^~ . In order to establish this, we derive an 
alternative approximation for Ao by starting from the identity (L0o,£o) — ^o(0o,£o) 
and making the approximation £o(u) ~ 1, which is valid outside a boundary layer 
around the absorbing boundary. Integration by parts now yields 

0o"(u*)V(u*)l v(w.)-0o(«*) ,,,m 

(00,1) (00,1) 

Moroever, from the initial condition (3.15) and the quasistationary approximation 
(3.23), we have 

(£o,p(O)) = (6)„ o (u_)~C o (So,0o). 
so that for £o( u ) ~ 1, 

(1,00) 

Equation (3.24) then shows that the first passage time density reduces to 

f(t) ~ A e- Aot (3.32) 

and(T)=/ oo t/(t)rft^l/A . 

3.2. WKB method and the quasi-stationary density. We now use the 

Wcntzel-Kramcrs-Brillouin (WKB) method [38, 52, 24, 47, 69] to compute an approx- 
imation e of 0o that does not satisfy the absorbing boundary condition. (Although 
such methods have been applied extensively to Fokker-Planck equations and master 
equations, the extension to CK equations of the form (3.9) is relatively recent). We 
thus seek an approximate solution of L<p e — of the WKB form 

e (w) = R(«) cxp f-^±\ , (3.33) 

where <!>(«.) is a scalar potential. Substituting into equation (3.21) gives 

(A + $'V)R = e(Vr)' + A R, (3.34) 

where ' denotes d/du. Introducing the asymptotic expansions R ~ R + eRi and 
3> = $o + e$i, and using the fact that Ao = 0{vT L /' i ), the leading order equation is 

AR = -$o VR o- (3-35) 

The diagonal components of V(w) are invertible almost everywhere for u > 0. Thus 
we can identify —$0 an d R-o as an eigcnpair of the eigenvalue problem 

M0 = /j,i/>, M = V _1 A. (3.36) 

Positivity of the probability density e requires positivity of the corresponding eigen- 
function 0. For fixed u, the matrix M has p as a right null vector and v as a left 
nullvector. These results follow from Ap = and v T V _1 A = 1 T A = 0. Thus one 
positive eigenfunction is O — P with corresponding eigenvalue fi = 0. However, such 
a solution is not admissible since $ = and $0 = constant. Since v n (u) for fixed 
< u changes sign as n increases from zero, theorem 3.1 of [58] ensures that there 
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exists one other positive eigenfunction, which we denote by xpi , whose corresponding 
eigenvalue /Xi varies with u in such a way that the corresponding WKB approximation 
is valid. Here we will construct such an eigenfunction explicitly. 

Using the explicit expressions for V and A, equations (3.10) and (3.11), the 
eigenvalue equation can be written in component form as 



n(—u + wn)ipr, 



(3.37) 



F(u)V„-i - (F(u) + n)ip n + (n+ 1)V„+1 
Trying a solution for ip\ of the form 

W»=7, (3.38) 

yields the following equation relating A and the corresponding eigenvalue /xi: 

■F(u) 



A 



-1 



■ + A — F(u) = /j,i(—u + wn). 



We now collect terms independent of n and linear in n to obtain the pair of equations 



1 



/'i 



F(u) 

A 



, A = F(u) — jiu. 



We hence deduce that 



and 



u 

-, Ml 

w 


1 

W 


wF(u) 


u 


(V'l)n = 


n 


\\w) 



Lb _L LU± \ LL ) , . 

A=-, /i 1 = ^-1 , (3.39) 



(3.40) 



where M is a normalization factor. Note that /ii(u) vanishes at the fixed points 
u_,u* of the mean-field equation (3.8) with /ii(u) > for < u < u_ and H\(u) < 
for m_ < u < u*. Moreover, comparing equation (3.5) with (3.40) establishes that 
ipi (u) = p(u) at the fixed points u. tl u±. In conclusion Ro = i/>i and the effective 
potential $o is given by 



(3.41) 



*o(«) = - / m(y)dy. 



The effective potential is defined up to an arbitrary constant, which has been fixed 
by setting $o( u -) = 0. 

Proceeding to the next order in the asymptotic expansion of equation (3.34), we 
have 

(A + fc'oVJRa = (VR )' - $;VR„. (3.42) 

Since Ro = ipi and $ = ~ Mi> it follows from the Fredholm alternative that 

T? T (V^l)' 



*i 



»7 T V^i 



(3.43) 
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where r\ is the left null vector of A — \i{W . Using equations (3.10) and (3.11), the 
components of r\ satisfy the explicit equation 

F(u)r) m+1 - (F(u) + m)r) m + mr\ m -\ = n\ \-u + mw]r) m . (3.44) 

Trying a solution of the form rj m — T m yields 

(F(u))T- (F^j + mj+mr 1 = m[-u + mw]. (3.45) 

r is then determined by canceling terms linear in to, which finally gives 



\w-b (u) 



Combining the various results, and defining 

k(u) = e*p(-J *i(y)dl/J, (3.47) 

gives to leading order in e, 



<f>,{«) -A7.-(i/)(-xi>( -^"Mj^c,,). (3.48) 



where ^2 n=0 ('4 > i)n — 1 and Af is a normalization factor, 



N 



k(u) exp 



$o(u) 



-l 



(3.49) 



/o 
The latter can be approximated using Laplace's method to give 



3.3. Perturbation analysis of the adjoint eigenfunction. Following Rcfs. 
[58, 41, 53], the adjoint eigenfunction £o can be approximated using singular per- 
turbation methods. Since Ao is exponentially small in e, equation (3.25) yields the 
leading order equation 

eV(u)^-+A T (u)t = Q, (3.51) 

supplemented by the absorbing boundary condition 

(€o)n(«*) = 0, n > k (3.52) 

A first attempt at obtaining an approximate solution that also satisfies the boundary 
conditions is to construct a boundary layer in a neighborhood of the unstable fixed 
point u* by performing the change of variables u = u* — ez and setting Q(z) = 
£o(w* — ez). Equation (3.51) then becomes 

V(u.)^+A T (u*)Q = 0. (3.53) 

dz 
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This inner solution has to be matched with the outer solution £o = 1, which means 
that 

lim QO) = 1. (3.54) 

z— >oo 

Recall from equation (3.36) that the (it-dependent) matrix M = V _1 A has eigen- 
values fij. Hence, introducing the similarity transform M = VMV" 1 and taking 
the transpose shows that M T = V~ 1 A T has the same eigenvalues. Denoting the 
corresponding eigenvectors by Cj we introduce the eigenfunction expansion 

oo 

Q(*) = col + X>Ci(«*)e« (u ' )z , (3.55) 

where 

^•(w)V(m)Cj(w) + A T (u)Cj{u) = 0. (3.56) 

In order that the solution remains bounded as z — > oo we require that Cj — if 
Pj{0) > 0. The boundary conditions (3.52) generate a system of linear equations 
for the coefficients Cj with codimension k. One of the unknowns in determined by 
matching the outer solution, which suggests that there arc k — 1 positive eigenvalues. 
The eigenvalues are ordered so that (J,j(0) > for j > k — 1. 

There is, however, one problem with the above eigenfunction expansion, namely, 
that ni(u*) — so that the zero eigenvalue is degenerate. Hence, the solution needs 
to include a secular term involving the generalized eigenvector £o; 

A T (u,)Co = -V(u*)l = -v(«*). (3.57) 

The Fredholm alternative ensures that £o exists, since p(u*) is the right null vector 
of A and p{u„) ■ v(u*) = 0, see equation (3.13). In component form with (Co)n — Qn, 

F(u*)(, n+ i + nCn-i - {F{u*) + n )Cn = u* — wn. (3.58) 

It is straightforward to show that this has the solution (up to an arbitrary constant 
that doesn't contribute to the principal eigenvalue) 

C« = wn. (3.59) 

The solution for Q(z) is now 

Q(z) = c l + ci« - zl) +Y j c J C J {u*> H{u ' )z - (3-60) 



J">2 



The presence of the secular term means that the solution is unbounded in the limit 
z — > oo, which means that the inner solution cannot be matched with the outer 
solution. One way to remedy this situation is to introduce an alternative scaling in 
the boundary layer of the form u = u* + e 1 ' 2 z, as detailed in Ref. [53]. Here we 
simply state the results of the analysis. The full inner solution takes the form 



&(«) 



u/ 



1/2 



i>2 
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The remaining coefficients c\, C2, • ■ • are determined by the boundary conditions (3.52), 
which reduce to 



nlW^^l + e^Coj^-gc^.)^! (3.62) 



for n> k. We thus find that 



c\ 



y 2 l^"*)l +0(^/2), c J =0(e 1 / 2 )forj>2 (3.63) 



3.4. Principal eigenvalue. It turns out that we only require the first coefficient 
c\ in order to evaluate the principal eigenvalue Ao using equation (3.29). This follows 
from the observation that Vrjj is an eigenfunction of the matrix A T V _1 , which are 
biorthogonal to the set of eigenvectors ipj of V _1 A. Since the WKB approximation 
4> e is proportional to i/>i, see equation (3.48), it follows that <p e is orthogonal to all 
eigenvectors C,, j ^ 1. Simplifying the denominator of equation (3.29) by using the 
outer solution £ ~ 1, we obtain 

6K) T vK)0eK) 

<^1> 



; Vi(m-)I / $oK) \ ,.,,... 

cifc(u*)B(u*)y — — — exp I — ), (3.64) 



with 



B(u t ) = Co" V ( u *)p( u *) = ^2 CnV n (u*)pn(u*) (3.65) 

n=0 

Substituting for ci and using the relation ^i(u) = $q(u), 

Ao ~ ifc(«*)B(tiO^(«_)|*o(«.)| exp f-^j , (3.66) 

Finally note that B(u*) can be evaluated using equations (3.5) and (3.59): 

oo 

£>(«*) = w > Pn(u*) \—u*n + wn ] 

= u> [— w*(n) + ui(n 2 )] . (3.67) 

Recall that p n (u) is given by a Poisson density with rate F{u), which implies that 
(n 2 ) = (n) + (n) 2 with (n) = ^(w*). Therefore, 

B(«*) = w [icF(u,) + F(!i»)(rf(«») -«*)]= f 2 %). (3.68) 

It is instructive to compare the effective potential $o obtained using the WKB 
approximation with the potential obtained using the FP equation (2.30) based on the 
QSS approximation. In the one-population case, equations (2.35) and (2.33) reduce 
to 

D = J2 «n(«)[«n(«) - V(u)] Pn (u), (3.69) 
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and 



J~] z n (u)A nm {u) = -[u m («) - V(x 



(3.70) 



Here v n and A are given by equations (3.10) and (3.11), p n is the Poisson distribution 
(3.7), and V(u) = —u + wF(u). At a fixed point, the equation for z n reduces to 
equation (3.57), and we find that z n — urn even away from fixed points. Substituting 
into equation (3.69) shows that 



D = (w 2 n[n - F(u)}) = w 2 F(u) = B{u) 



(3.71) 



with B given by equation (3.68). The steady-state solution of the FP equation (2.30) 
takes the form C{u) ~ exp~*°(")/ c with stochastic potential 



I D{y) y J w*F(y) 



(3.72) 



Note that <I> differs from the potential $o, equation (3.41), obtained using the more 
accurate WKB method. Equations (3.39) and (3.41) show that the latter has the 
integral form 



$o(u) = - 



u l 



■>F(y) 
V 



- 1 



dy. 



(3.73) 



Thus, there will be exponentially large differences between the steady-states for small 
e. However, it gives the same Gaussian-like behavior close to a fixed point u*, that is, 



du 



du 



0, 



a 2 $ 



du 2 



d 2 <P 



du 2 



1 - wF'(u) 



U — U^ 

(3.74) 



0.40 




Fig. 3.2. Comparison of the double-well potentials $o(i*) o,nd 3>o(") obtained using the quasis- 
tationary approximation and the QSS diffusion approximation, respectively. Parameter values are 
chosen so that deterministic network is bistable: Fq = 2, ■y = 4, k = 1, and w = 1.15. 
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3.5. Results. In Fig. 3.2, we plot the potential function $o of equation (3.73), 
which is obtained using the quasistationary approximation in a parameter regime for 
which the underlying deterministic network is bistable. We also plot the correspond- 
ing potential function <£>o of equation (3.72), under the QSS diffusion approximation. 
The differences between the two lead to exponentially large differences in estimates 
for the mean exit times when e is small. The mean exit time from the left and right 
well is shown in Fig. 3.3. Solid curves show the analytical approximation T ~ 1/Ao, 
where Ao is given by (3.64), as a function of 1/e. For comparison, the mean exit time 
computed from averaged Monte-Carlo simulations of the full stochastic system are 
shown as symbols. From (3.64), we expect the log of the mean exit time to be an 
asymptotically-linear function of 1/e, and this is confirmed by Monte-Carlo simula- 
tions. The slope is determined by the depth of the potential well, and the vertical shift 
is determined by the prefactor. Also shown is the corresponding MFPT calculated 
using the QSS diffusion approximation (dashed curves), which is typically several 
orders of magnitude out, and validates the relative accuracy of the quasistationary 
approximation. 

4. Metastable states in a two population model. The same basic analyti- 
cal steps as §3 can also be used to study multipopulation models (M > 1). One now 
has M piecewise-deterministic variables U a and M discrete stochastic variables N a 
evolving according to the CK equation (2.19). The mean first passage time is again 
determined by the principal eigenvalue Ao of the corresponding linear operator. As 
in the one-population model, Ao can be expressed in terms of inner products involv- 
ing a quasi stationary density (fi e , obtained using a multidimensional WKB method, 
and the principal eigenvector £ of t ne adjoint linear operator, calculated using sin- 
gular perturbation theory. One of the major differences between the one-population 
model and multi-dimensional versions is that the latter exhibit much richer dynamics 
in the mean-field limit, including oscillatory solutions. For example, consider a two- 
population model (M = 2) consisting of an excitatory population interacting with 
an inhibitory population as shown in Fig. 4.1. This is one of the simplest determin- 
istic networks known to generate limit cycle oscillations at the population level [5], 
and figures as a basic module in many population models. For example, studies of 
stimulus induced oscillations and synchrony in primary visual cortex often take the 
basic oscillatory unit to be an E-I network operating in a limit cycle regime [70, 36]. 



io 5 



h: 



io 1 



left well 



1/e 



right well 




FlG. 3.3. Mean exit time from the left and right well calculated using the quasistationary 
approximation (solid line) and the QSS diffusion approximation (dashed line). The open circles 
represent data points obtained by numerically solving the corresponding jump velocity Markov process 
using the Gillespie algorithm. Parameter values are the same as in Fig. 3.2. 
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Here the E-I network represents a cortical column, which can synchronize with other 
cortical columns either via long-range synaptic coupling or via a common external 
drive. In this paper, we will focus on parameter regimes where the two-population 
model exhibits bistability. 




Fig. 4.1. Two-population E-I network with both intrasynaptic connections Wee,Wii and 
intersynaptic connections Wie,Wei- There could also be external inputs hE,hj, which can 
be incorporated into the rate functions of the two-populations by shifting the firing threshold 

Let u\ = x and U2 — y denote the synaptic variables of the excitatory and 
inhibitory networks, respectively, and denote the corresponding spiking variables by 
n x and n y . The CK equation (2.19) can be written out fully as 

dp d(vp) d(vp) 

dt dx dy 

1 



(4.1) 



+ - [w-(n x + l)p(x,y,n x + l,n y ,t) +u>-(n y + l)p(x,y,n x ,n y + l,i)] 
+ - [^+(x)p(x ) y,n x - l,n y ,t) + u + (y)p(x,y,n x ,n y - l,i)] 



where 



and 



[u-(n x ) +u-{n y ) + u+(x) + uj + (y)]p{x,y,n x ,n y ,t), 



v(x, n x , n y ) = -x + [w E En x - w E in y ] 



v(y, n x , n y ) =-y+ [w IE n x - w n i 



uj+{x) — F(x), w_(n) 



ui ■ 



(4.2) 
(4.3) 

(4.4) 



Thus the synaptic coupling between populations occurs via the drift terms v, v. As 
in the case of the one-population model, we expect the finite-time behavior to be 
characterized by small Gaussian fluctuations about the stable steady-state of the 
corresponding pure birth-death process. We now show that in the limit N — > oo 
and At —> with NAt = 1 and e fixed, the steady-state distribution reduces to a 
multivariate Poisson process. First, introduce the generating function (for fixed (x, y)) 



G(r,s)= Y^ Yl r nx s n «p(x,y,n x ,n v ). 



(4.5) 
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Setting all derivatives in equation (4.1) to zero, multiplying both sides by r" x and r"» 
and summing over n x , n y gives the quasilinear equation 



= (1-r) 



dG . ,dG 
/') — +(l-s) 



+ [(r-l)w + (s) + (a-l)w+(y)]G. 



dr ds 

This can be solved using the method of characteristics to give 

G(r, s) = exp ([r - l]w+(x) + [s - l]w+(y)] , 

which is the generating function for the steady-state Poisson distribution 



p(x,y,n x ,n y ) 



I 



n y ! 



Since (n x ) = w+(a;), (n y ) = w_(y), it immediately follows that in the limit e 
obtain the standard voltage-based mean-field equations for an E-I system: 

dx 

~dt 

dy 
dt 



■y - 



- w E eF(x) 
w IE F(x) 



-w EI F(y), 
w n F(y) 



(4.6) 
(4.7) 

(4.8) 
> 0, we 

(4.9) 
(4.10) 



It is well known that the dynamical system (4.9) exhibits multistability and limit 
cycle oscillations [5]. We will assume that it is operating in a parameter regime for 
which there is bistability as illustrated in Fig. 4.2. For example, if wee — Wei = 
Wie — wii — w then x = y is an invariant manifold on which x evolves according to 
the one-population equation (3.8). Varying the threshold n then leads to a pitchfork 
bifurcation and the emergence of bistability. 
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Fig. 4.2. Deterministic limit of the two population model, showing bistability. Red curves 
show the x-nullclines, and blue curve show the y-nullcline. The red nullcline through the 
saddle is its stable manifold and acts as the separatrix T between the two stable fixed points. 
Two deterministic trajectories are shown (black curves), starting from either side of the 
unstable saddle and ending at a stable fixed point. Parameter values are Fq — 1, 7 = 3, 
K = 2, wee = 5, Wei = 1, Wie = 9, and wn = 6. 

In order to analyze the effects of fluctuations for < e < 1, we rewrite equation 
(4.1) in a more compact form by introducing some tensor notation. First, we introduce 
the probability 1-tensor p(x,y,t) with components 



Pr, 



,{x,y,i) =p(x,y,n x ,n v ,t) 



(4.11) 
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the diagonal drift 2-tensors V(x), V(y) with diagonal components 

Vn x ,n y ;n x ,n y (%) =v(x,n x ,n y ), Vn x ,n y ;n x ,n y (l/) = v(y, Tl x , fly), (4-12) 

and the transition 2-tensor A(x, y) with non-zero components 

An x ,n y ;n x -l,n y (x,y) = F(x), A^^.^^^X.y) = F(j/), (4.13) 

An x ,n y ;n x +i,n y {x,y)=n x + l, A nir . ny . nx ^ y+ i (x, y) = n y + 1, (4-14) 

and 

A^.n^.n.fos/) = -[F(x) + J F(y) + n a; + nj / ]. (4.15) 

Second, we rewrite the CK equation as 

|^-A(Vop)-|-(Vop) + iAop, (4.16) 

where 

[■"- \?\n x ,n y / ^na; ,n y ;m x ,m y Vm x ,m y •> V '/ 

m x ,m y 

etc. The tensor A satisfies the null constraints (cf. equation (3.12) 

l T oA = 0, Aop = 0, (4.18) 

with l nxi „ y = 1 for n x ,n y and pn x , ny (x,y) = p(x,y,n Xl n y ) with the density p given 
by the Poisson distribution (4.8). In tensor notation, the mean-field equations (4.9) 
can be written as 

dx dv ~ , 

— =Vop, J-Vop. (4.19) 

Finally, note that as in the one-population model, we can restict the domain of the 
stochastic dynamics in the (x, ?/)-plane. In order to show this, multiply both sides of 
equations (4.2) and (4.3) by wu and wei respectively, and add the resulting equations. 
Setting 

x = [wjix — w e iy] / det[w] (4.20) 

with det[w] = weeWu — WeiWie, we have the transformed drift term 

v{x,n x ) = -x + n x . (4-21) 

Similarly, multiplying both sides of equations (4.2) and (4.3) by wie and Wee respec- 
tively, and adding the resulting equations yields 

v(y,n v ) = -y + riy. (4.22) 

with 

y = [wiex - weev]/ dek[w]. (4.23) 

It follows that the dynamics can be restricted to the domain x > 0, y > 0. 
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4.1. Quasi-stationary approximation. In order to investigate rare transitions 
between the metastable states shown in Fig. 4.2, we introduce an absorbing boundary 
along the separatrix T separating the two states: 

p(x,y,t)=0, (x,y)£T (4.24) 

for all components (n xi n y ) for which 

(v(x, n x , n y ),v(y, n x ,n y )) • s < 0, (4-25) 

where s is the unit normal of T pointing into the domain T> of the initial metastable 
state, which we take to be (x_,y_). Following identical arguments to the one- 
population model, we can expand the probability density as 

DO 

p(x,y,t) =^c,0,(a;,y)e- A ^ (4.26) 

3=0 



with (Xj,(pj) determined from the eigenvalue equation 
together with the boundary conditions 



Lfy = — (V o <Pj) + — (V o 0,) - - A o fy = Xj4>j , (4.27) 



(<pj)(x,y,n x ,n y ) = 0, for all(.x,y) e T, a,nd(n x ,n y ) (4.28) 

for which equation (4.25) is satisfied. The principal eigenvalue Ao again determines the 
first passage time density according to f(t) ~ A c~ A °'. Moreover, A can be approx- 
imated using a spectral projection method that makes use of the adjoint eigenvalue 
equation 



!(*)-?„ » (W -i. 



L "d 5-Vo-ft)-V» F ft|--Ao{ i =A^, (4.29) 



with (4>i,£j) = Si j and the inner product defined for 1-tensors according to 

(f,g)= / f{x,y) T og{x,y)dxdy, 
Jv 

= y~] / f(x,y,n x ,n y )g(x,y,n x ,n y )dxdy. (4.30) 

Now suppose that we replace </> by the quasi-stationary density (p e , for which L<fi e = 
without satisfying the absorbing boundary conditions. Application of the divergence 
theorem shows that 

(<p e , L*£ ) = (L<p € , £ ) + J ($ o V o e , $ v ^ . zds. (4.31) 

It follows that 

^°Vo4,5 T oVo e ) • rids 

. (4.32) 



J v ^o(f> e dA 
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4.2. WKB method and the quasi-stationary density. Following along sim- 
ilar lines to the one-population model, we approximate the quasi-stationary density 
4> e (x,y) of the CK equation (4.16) using the WKB method. That is, we seek an 
approximate solution of L<p e — of the WKB form 

ej> e (u) = {no(x,y) + eR 1 (x,y)]e X p(- * o{x ' y) + e * l{x ' y) ), (4.33) 



Here R and Hi are 1-tensors and $01*^1 are scalars. Substituting into equation 
(4.16) and collecting leading-order terms in e gives 

[A + V X V + T y V] o R = 0, (4.34) 

where 

Given the explicit form of the diagonal tensors V, V, see equation (4.12), P{V + V2V 
for (x, y) G T> has at least two components of opposite sign. This is a necessary 
condition for the existence of a non-trivial positive solution for Ro in the domain T> 
according to Theorem 3.1 of [58] . 
We now make the ansatz 

[Rok,„ y = "fr • -T- (4-36) 

n x \ n y \ 

Substituting into equation (4.34) and using the explicit expressions for A, V and V, 
we find that 



F(x) 1 



F(y) 1 



Ay 



i y + A x +A y - F(x) - F{y) 



A x 

= -V x \-x + w EE n x - w E in y ] - V v [-y + w IE n x - Wnn y ]. (4.37) 



The variables V x and V y can be determined by cancelling terms in n x and n y . This 
yields the pair of simultaneous equations 

Fix) F(v) 

-j- 1 - 1 = -[weeV x + wiEVy}, -^ - 1 = w EI V x + wiiVy. (4.38) 

x y 

Substituting back into equation (4.37) gives 

xV x + yV v = A x + A y - F(x) - F(y). (4.39) 

Solving for A x , A y in terms of V x and V y , equation (4.39) can be rewritten as 

U{x,y,V x ,Vy) = -xV x - yV y - F(x) - F (y) + A x {x,V x ,V y ) + A y {y,V x ,V y ) = 0, 

(4.40) 
where 



F(x) A F(y) 

1 - w EE V x - w IE V y ' y 1 + weiFx + wuPy 



K = , - ^' - „ , A, = V ■ ... -r, ( 4 - 41 ) 



Mathematically speaking, equation (4.40) is identical in form to a stationary Hamilton 
Jacobi equation for a classical particle moving in the domain T>, with % identified as 
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the Hamiltonian. A trajectory of the particle is given by the solution of Hamilton's 
equations 



dx 


dH dy dH 


dt 


' dV x ' dt ~ dV y ' 


dV x 


dH dV y dH 


dt 


dx ' dt dy 



(4.42) 

Here t is treated as a parameterization of trajectories rather than as a real time 
variable. Given a solution curve (x(t),y(t)), known as a ray, the potential $0 can be 
determined along the ray by solving the equation 

d$ <9$n dx <9$ dy _ dx _ dy , 

dt dx dt dy dt dt y dt 

Thus, $0 can be identified as the action along a zero energy trajectory. One can then 
numerically solve for $0 by considering Cauchy data in a neighborhood of the stable 
fixed point (x_ , y_ ) [53] . We find that Hamilton's equations take the explicit form 

dx weeV x + w IE Vy „, . 

— = -x + weeF(x) - WEiF(y) + T ^ttoWeeF(x) 

dt 1 - weeV x - w IE Vy\ 2 



WElVx + WllVy 

[weiV x + w H V v + 1] 



;W E iF(y) (4.44) 



dy , . , . weeVx + wi E V v 

— = -y + w IE F(x) - WiiF{y) + — -w IE F{x) 

at [1 - weeVx - wi E V y \ l 

WElVx+WuVy 

jWnF(y) (4.45) 



[wEiVx + wnVy + l]' 



d ^=Vx- WeeV * + WieV v F > {x) (4.46) 

dt x l-WE E Vx-w IE Vy K ' v ' 



®\_ v WElVx+WuVy l() , . 

dt -' y+ weeVx + w IE V y + r {V > (4 ' 47 > 



Note that we recover the mean-field equations along the manifold V x = V y = with 
Ax = F(x),A y = F(y). 

It remains to specify Cauchy data for the effective Hamiltonian system. At the 
stable fixed point, the value of each variable is known with V x — V y = and (x, y) = 
(x_, y_). However, data at a single point is not sufficient to generate a family of rays. 
Therefore, as is well known in the application of WKB methods [52, 47, 69], it is 
necessary to specify data on a small ellipse surrounding the fixed point. Thus, Taylor 
expanding $0 around the fixed point yields, to leading order, the quadratic form 

^(x,y)^ 1 -z T Zz, z=( X ~ X - V (4.48) 



where Z is the Hessian matrix 



2 ' V y-y. 



3 2 $o d 2 $o 

dx 2 dxdy 

d 2 $o d 2 $ 

dydx dy 2 



(4.49) 
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and we have chosen $o(^-, 2/-) = 0. Cauchy data are specified on the ^-parameterized 
ellipse 

±z T (8)Zz(6) = X , (4.50) 

for a suitably chosen parameter x that is small enough to generate accurate numerical 
results, but large enough so that the ellipse can generate trajectories that cover the 
whole domain T>. On the elliptical Cauchy curve, the initial values of V x and V y are 

v vfi (0) )- Z {yo(e)-y- )■ (451) 

It can be shown that the Hessian matrix satisfies the alebraic Riccati equation [47] 

ZBZ + ZC + C T Z = 0, (4.52) 

where 

/ §Ph d 2 u \ ( d 2 u d 2 u \ 

R | dV * dV^dVy \ r _ I dVJix dV f dy \ (A con 

\ 97^97^ 9PJ / \ dV y dx dV y dy ) 

are evaluated at V x — V y — and (x, y) = (x_,y_). In order to numerically solve 
the Ricatti equation it is convenient to transform into a linear problem by making the 
substitution Q = Z _1 : 

B + CQ + QC T = 0. (4.54) 

Proceeding to the next order in the WKB solution of equation (4.16), we find 
that 

Since the 2-tensor M = A + V X "V + VyV has the unique right null 1-tensor Ro, it's 
left null-space is also one-dimensional spanned by rj, say. The Fredholm alternative 
theorem then requires that 



T 
r\ o 



d[V°R ] + g[v°R ] _ /a*i v + 5*1 v . o R 

dx dy \ dx dy 



(4.56) 



Using the fact that (r) T o V o H )dy/dt — (r/ T o V o R Q )dx/dt along trajectories of 
the Hamiltonian system, we can rewrite the above equation as (cf. equation (3.43)) 

d<$> 1 _d<5> l dx d^xdy _ dx/dt T /<9[VoR ] <9[VoR ]\ 

~dT = ~fa~~di + ~dy"di ~ rf oV oRq 77 ° ^ dx + dy )' ( ' ' 

As shown in appendix A of [53] , an equation of this form can be numerically integrated 
along the trajectories of the underlying Hamiltonian system. 

However, rj may be solved explicitly by substituting the ansatz 

Vn x ,n y =T^-T^ (4.58) 
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into the equation rj o M = 0, and using the explicit expressions for A, V, V. One 
finds that 



F(x)[T x -l}+F(y)[T y -l] 



1 



- 1 



1 



- 1 



= -V x \-X + w EE n x - WEin y ] - P y [-y + w IE n x - w u n y ]. 
Cancelling the terms in n x and n y yields 

1 



r, 
l 



- 1 = -[w E eV x + wi E Vy], 



- 1 = WElVx + WllVy 



Comparison with equation (4.38) shows that 

A, „ A,, 



r„ 



so that 



>)„ 



F(,-V '» F(y) 
A, \" f Ay 



(4.59) 

(4.60) 
(4.61) 

(4.62) 
(4.63) 



F(x)J \F(y) 

In summary, the quasi-stationary approximation takes the form 

<t>e(x,y) ~ N e -* l( - x 'rte-* o( - x ' y V e Ro. (4.64) 

The normalization factor M can be approximated using Laplace's method to give 

-l 



N = 



exp 



v 



$o{x,y) 



*i(*,y) 



y/det(Z(x-,y-)) 
2ttc 



where Z is the Hessian matrix 



3 2 $o d 2 $ \ 

dx 2 dxdy \ 

d 2 $o d 2 <S> I ' 

dydx dy 2 / 



(4.65) 



(4.66) 



and we have chosen <&o(x-, y_) = = <&i(a;_,y_). 

4.3. Perturbation analysis of the adjoint eigenfunction. Since Ao is expo- 
nentially small, the leading-order equation for the adjoint 1-tcnsor to is 



dx dy 



o £o + A T o £o = 0, 



(4.67) 



supplemented by the absorbing boundary conditions (with (£o)n x ,n — £,n x .n ) 

Zn x ,n y (x,y) = 0, (x,y)eT (4.68) 

for all (n x ,n y ) such that 

(v(x, n x , n y ),v(y, n x ,n y )) • s > 0. (4.69) 
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Following along similar lines to [53], wc introduce a new coordinate system r = 
r(x,y),a — a(x,y) such that r parameterises the separatrix (x,y) G T and a is a 
local coordinate along the normal s of V. We scale a so that (dx/da, dy/da) = s at 
a = 0. Equation (4.67) becomes 



V r 



d_ 

~frr 



V„ 



d 



o £ + A T o £ = 0, 



where 



dx 



Or, 

dy 



V. 



ox 



dy 



V. 



(4.70) 



(4.71) 



and all terms are rewritten as functions of r, a. Thus, A(er, r) = A(a;(<7, r),y(a, r)) 
etc. As a first attempt at obtaining an approximation for £oi we introduce a boundary 
layer around T by setting a — ez and Q(z,r) = ^o( e - z 7 T )- To leading-order, equation 
(4.70) becomes 



V„(0,r) 






Q(z, r) + A T (0, r) o Q(z, r) = 0. 



(4.72) 



The inner solution has to be matched with the outer solution ^ = 1, which means 

lim Q(x,t) = 1. (4.73) 

We now introduce the eigenfunction expansion (cf. equation (3.55)) 

Q(z,r) = C0 (T)l + ]T Cj (T)O(0,T)e^ ^, (4.74) 

i>i 

where 1 has a zero eigenvalue, and 

Mj (a,T)V ff (a,T)oC,( ( T,T) + A T (a,T)oC J (a,T) = 0, j^0. (4.75) 

In order that the solution remain bounded as z — » oo and fixed r, we require that 
Cj(r) = if fij(0,r) > 0. Suppose that the boundary conditions (4.68) for fixed r 
generate a system of linear equations for the unknown coefficients Cj (r) of codimension 
k. One of the coefficients is determined by matching the outer solution, which suggests 
that there are k — 1 positive eigenvalues for each r. The eigenvalues are ordered so 
that for each r, /Uj(0, r) > for j > k — 2. 

Analogous to the one-population model, an additional eigenvalue \x\ , say, vanishes 
at the saddle point (0,r*) on the separatrix. In order to shows this, suppose that 



t'i 



1 T o [V X V + V y V] o Ci = d^p 
l T oV„o Ci da ' 



(4.76) 



where the last expression follows from equations (4.35) and (4.71). Substitution into 
equation (4.75) for j = 1 then gives 



[A T + V X V + V y V] o Ci = 0, 



(4.77) 



which has the unique solution £i = rj, the adjoint of Ro Since V x and V y vanish at 
(0, r*) and Vg- o r\ ^ 0, it follows that /tti(0, r*) = 0. Hence, the solution at r* has to 
include a secular term involving the generalized eigenvector £qj where 



A T (0, n ) o Co = -V ff (0, r, ) o l = -s • (V, V) o l. 



(4.78) 
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The Frcdholm alternative theorem ensures that a solution exists, since p(0, r*) is the 
left null 1-tcnsor of A t (0,t*) and p(0,r*) o V(0,r*) = p(0,r*) o V(0,r*) = 0. More 
explictly, setting [Co]n x ,n„ = C{n x ,n y ), we have 

F{x)Q{n x + 1, rij,) + F(y)((n x ,n y + 1) + n x ((n x - 1, n„) + n y C,{n x ,n y - 1) 

~[F(x) + F(j/) +71^+ %]((«-£, %) 

= s-Jx - w EE n x + w EI n y ] + s y [y - w IE n x + w n n y ]. (4.79) 

This has a solution of the form £„ x , n = yt^ria; + A y n y , with the coefficients *4a;,„4y 
determined by canceling linear terms in n x ,n y . Thus 

Cn x ,n„ = [w EE s x + w IE s y ]n x - [w EI s x + w H s y ]n y . (4.80) 

Given £o, equation (4.74) becomes 

Q(z,n) = c (n)l + Cl (n)(Co - zl) + Y, c ^KMr*)^ ( ^ )z , (4.81) 

The presence of the secular term implies that the solution is unbounded so it has 
to be eliminated using a modified stretch variable a = yfez [54, 53]. As in the one- 
population case, we find that 



Cl(T .) „ _^e|cW0,T.)| (4.82) 

4.4. Principal eigenvalue. We now return to the expression for the principal 
eigenvalue Ao given by equation (4.32). Simplifying the denominator by using the 
outer solution £o ~ 1 and using the WKB approximation of <p e , equation (4.64), gives 

A = N f e-* 1 ^'»)e-* o(a, ' , ' )/e (^ ° V o R , jJoVoRj). Sds. (4.83) 

Changing to the (a, r) coordinate system and evaluating the line integral by applying 
Laplace's method around the saddle point (0, r*) then gives 

A ~MB(T*) Cl (T*)e-^ ^e-*°^/*[ exp f-^ rT $o(0,r,)(r-T,) 2 j rfr, 
^i?(n) Cl (n)e-^^)e-^^)/%/ n 2 " " AAMZQc-.y-)) 

k ;n ; ya Tr $ (o,T*) 2^ e 

~ 1 g(- r Jc-*i^')c-*''^ T '^%/ ^' T ' T$o(0 ' T * )det(Z(a: -' y - )) , (4.84) 

7T * y <9 Tr $0(0,1-*) 

where we have used equations (4.65), (4.81), (4.82), and 

B(n) = U% o V o R o (0, t.), (o T| >Vo Ro(0, r.)) • n. (4.85) 



4.5. Results. The rays (x(t),y(t)) (i.e., solutions to the Hamilton's equations 
(4.44) in the (x, j/) plane) have an important physical meaning. The trajectory of 
the ray is the most likely trajectory or path leading away from a stable fixed point 
[24]. Under this interpretation, one can describe the stochastic dynamics using the 
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Fig. 4.3. (a) Characteristic paths of maximum likelihood for the 2D model. Rays originating 
from the left (right) stable fixed point are shown in orange (cyan), with the ray connecting to the 
saddle shown in red (blue). The grey curve is the separatrix Y . Level curves of constant <&(£, y) are 
shown as black dots. Each ray has four dots for different values of &(x,y). Rays originating from 
the left fixed point have dots at $ = 0.1, 0.2, <£>* + 0.01, $« + 0.02, and rays originating from the right 
fixed point have dots at 5> = 0.19,0.23,0.28,0.30, where 5>* = &(x,,y t ) = 0.28. All rays terminate 
at $ = <t>*+0.02. (b) Sample trajectories of the two-population velocity jump Markov process, whose 
associated probability density evolves according to (4-1), are computed using the Gillespie algorithm 
with e = 0.05 and NAt = 1. (The maximum likelihood paths are independent of epsilon). Other 
parameter values are the same as in Fig. 4.2. 



metastable dynamical trajectories (rays) along with deterministic trajectories. The 
rays (x(t),y(t)) shown in Fig. 4.3 are obtained by integrating the characteristic equa- 
tions (4.44) and (4.45). These trajectories are only valid in one direction: away from 
the stable fixed points. The most likely trajectory leading toward stable fixed points 
are given by deterministic dynamics (see Fig. 4.2). For parameter values considered 
in Fig. 4.3, rays originating from each stable fixed point cover separate regions, so 
that most likely paths between points in each region are connected by deterministic 
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trajectories starting at the boundary between the two regions. Note that this bound- 
ary is not the separatrix (grey curve). For example, a trajectory initially at the left 
fixed point which crosses the separatrix at the saddle would most likely follow a ray 
toward the saddle and then follow a deterministic trajectory to the right fixed point. 
If a trajectory crosses the separatrix away from the saddle, it is most likely to cross 
the separatrix above the saddle when starting from the left fixed point and below the 
saddle when starting from the right fixed point (see Fig. 4.4). At first glance, this 
suggests that if the trajectory starts at the left fixed point, say, it is more likely to 
cross above the saddle, continue along a deterministic trajectory to the right fixed 
point, and then cross the separatrix below the saddle than it is to directly cross below 
the saddle. This is counter to intuition because it would seem more likely for a single 
rare, metastable crossing event to lead to a point near the separatrix than two rare 
events occurring in sequence. However, as shown in [47], rays can also originate from 
the saddle point that cross the separatrix in the direction oposite those originating 
at the stable fixed points. In Fig. 4.5, the probability density function for the y 



y 




Fig. 4.4. Maximum- liklihood trajectories crossing the separatrix. 



coordinate of the point on the separatrix reached by an exit trajectory is shown for 
each well (square symbols show the histogram for exit from the left well and likewise, 
'o' symbols for the right well). Each density function is peaked away from the saddle 
point, showing a phenomena known as saddle point avoidance [47, 69]. As e — > 0, the 
two peaks merge at the saddle point. Although we expect the saddle point to be the 
most likely exit point since it the point on the separatrix where the potential $ takes 
its minimum value, our results show that this is not necessarily true. 

Even though the most likely exit point is shifted from the saddle, the value of 
potential at the saddle point still dominates the mean first exit time. In Fig. 4.6, the 
mean exit time from each of the 2D potential wells (see Fig. 4.3) is shown. Solid lines 
show the analytical approximation T ~ 1/Ao, where Ao is given by (4.84), and symbol 
show averaged Monte-Carlo simulations. As in Fig. 3.3, the slope T on a log scale as 
a function of 1/e is determined by <I> evaluated at the saddle point. 

5. Discussion. In this paper we developed a generalization of the neural master 
equation [17, 7, 18], based on a velocity jump Markov process that couples synaptic 
and spiking dynamics at the population level. There were two distinct time-scales in 
the model, corresponding to the relaxation times r and r a of the synaptic and spiking 
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Fig. 4.5. The probability density for the exit point (y coordinate) where the separatrix is crossed 
by an exiting trajectory. Results are obtained by 10 2 Monte-Carlo simulation with the same param- 
eters as used in Fig. J^.2, with e = 0.08. The square symbols show trajectories from the left well, and 
'o' symbols show trajectories from the right well. 
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Fig. 4.6. Mean exit time from the left and right well. Parameter values are the same as in 
Fig. 4-Z. Solid lines show the analytical approximation T ~ 1/Ao, where Ao is given by (4.84), and 
symbol show 80 averaged Monte-Carlo simulations. 



dynamics, respectively. In the limit r — > 0, we recovered the neural master equation 
operating in a Poisson-like regime, whereas in the limit t„ -> we obtained determistic 
mean field equations for the synaptic currents. Hence, one additional feature of our 
model is that it provides a prescription for constructing a stochastic population model 
that reduces to a current-based model, rather-than an activity-based model, in the 
mean-field limit. 

We focused on the particular problem of escape from a metastable state, for which 
standard diffusion-like approximations break down. We showed how WKB methods 
and singular perurbation theory could be adapted to solve the escape problem for 
a velocity jump Markov process, extending recent studies of stochastic ion channels. 
For concreteness, we assumed that the network operated in the regime T a /r = 6 « 1, 
which meant that transitions between different discrete states of population spiking 
activity were relatively fast. In this regime, the thermodynamic limit N — > oo was 
not a mean- field limit, rather it simplified the analysis since the quasi-steady-state 
density was Poisson. It would be interesting to consider other parameter regimes in 
subsequent work. First, we could model the discrete Markov process describing the 
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spiking dynamics using the Bressloff version of the master equation [7] . There would 
then be two small parameters in the model, namely e and 1/iV, so one would need to 
investigate the interplay between the system size expansion for large but finite N and 
the quasi-steady-state approximation for small e. Another possible scenario (though 
less plausible physiologically speaking) would be fast synaptic dynamics with t <C T a . 
In this case, mean-field equations are obtained in the thermodynamic limit. Finally, 
it would be interesting to extend our methods to analyze the effects of noise when the 
underlying deterministic sytstem exhibite more complicated dynamics such as limit 
cycle oscillations. As we commented in the main text, the two-population model of 
excitatory and inhibitory neurons is a canonical circuit for generating population-level 
oscillations. 

Finally, it is important to emphasize that the neural master equation and its gen- 
eralizations are phenomenological models of stochastic neuronal population dynamics. 
Although one can give a heuristic derivation of these models [9] , there is currently no 
sytematic procedure for deriving them from physiologically-based microscopic models, 
except in a few special cases. Nevertheless, stochastic hybrid models are emerging in 
various applications within neuroscience, so that the analytical techniques presented 
in this paper are likely to be of increasing importance. 
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